

library(matchMulti)
library(dplyr)
library(ggplot2)

# Load Data
rm(list=ls())
data(catholic_schools)

# Desc Stats
describe_data_counts( catholic_schools, "school", "sector" )

catholic_schools$sectorf <- factor(catholic_schools$sector, label=c("Public", "Catholic"))

#Some Box Plots to Examine Common Support
ggplot(catholic_schools, aes(x=sectorf, y=acad)) + 
		geom_boxplot() + guides(fill=FALSE) + xlab("") + 
		ylab("Percentage of Students on Academic Track") +
		scale_fill_brewer(palette = "Paired") +
       theme(text = element_text(size=20)) 

ggplot(catholic_schools, aes(x=sectorf, y=discrm)) + 
		geom_boxplot() + guides(fill=FALSE) + xlab("") + 
		ylab("Disciplinary Climate") +
		scale_fill_brewer(palette = "Set2") +
        theme(text = element_text(size=20)) 

ggplot(catholic_schools, aes(x=sectorf, y=female_mean)) + 
		geom_boxplot() + guides(fill=FALSE) + xlab("") + 
		ylab("Percentage of Female Students") +
		scale_fill_brewer(palette = "Set3") +
        theme(text = element_text(size=20))
		 
# Secondary Plots 
ggplot(catholic_schools, aes(x=sectorf, y=size, fill=sectorf)) + geom_boxplot() +
  guides(fill=FALSE) + xlab("") + ylab("School Enrollment") 
  
ggplot(catholic_schools, aes(x=sectorf, y=minority_mean, fill=sectorf)) + geom_boxplot() +
  guides(fill=FALSE) + xlab("") + ylab("Percentage Minority") 


data <- catholic_schools %>% filter(female_mean>.30, female_mean<.75)

# Overlap Diagnostic
school.data <- data %>% group_by(school) %>% 
   summarize(across(c('minority_mean', 'female_mean', 'size', 'acad', 'discrm', 'ses_mean', 'sector'), ~mean(.x))) %>% as.data.frame()
                
pscore <- glm(sector ~ minority_mean + female_mean + size + acad + discrm + ses_mean, data=school.data)
school.data$lin.ps <- pscore$linear.predictor

trt <- school.data %>% filter(sector==1)
ctrl <- school.data %>% filter(sector==0)

diff <- outer(trt$lin.ps, ctrl$lin.ps, FUN = "-")
diff <- abs(diff)

test <- diff < .10
zeta.t <- apply(test, 1, sum)
sum(zeta.t > 1)/32

### 50% of Catholic Schools have close comparisons in the control group.

test <- diff < .05
zeta.t <- apply(test, 1, sum)
sum(zeta.t > 1)/32


